﻿using LatoolNet;

namespace NumericalSolution {
  class ConvectionDiffusion1DCentral : ConvectionDiffusion1D {
    public override Matrix Solve() {

      int bandWidth = 3;
      Matrix mat = new Matrix(p_numberOfControlVolumes,
                              p_numberOfControlVolumes, bandWidth);
      Matrix vec = new Matrix(p_numberOfControlVolumes, 1);

      double phi_0 = 1.0;
      double phi_L = 0.0;

      mat[0, 0] = (p_rho * p_u) / 2.0 + 3.0 * p_gamma / p_deltax;
      mat[0, 1] = (p_rho * p_u) / 2.0 - p_gamma / p_deltax;
      vec[0, 0] = (p_rho * p_u) * phi_0 + 2.0 * p_gamma * phi_0 / p_deltax;

      for (int i = 1; i < p_numberOfControlVolumes - 1; i++) {
        mat[i, i - 1] = - (p_rho * p_u) / 2.0 - p_gamma/ p_deltax;
        mat[i, i] = 2.0 * p_gamma / p_deltax;
        mat[i, i + 1] = (p_rho * p_u) / 2.0 - p_gamma / p_deltax;
        vec[i, 0] = 0;
      }

      mat[p_numberOfControlVolumes - 1, p_numberOfControlVolumes - 1]
        = - (p_rho * p_u) / 2.0 + 3.0 * p_gamma / p_deltax;
      mat[p_numberOfControlVolumes - 1, p_numberOfControlVolumes - 2]
        = - (p_rho * p_u) / 2.0 - p_gamma / p_deltax;
      vec[p_numberOfControlVolumes - 1, 0]
        = - (p_rho * p_u) * phi_L + 2.0 * p_gamma * phi_L / p_deltax;

      //LUFactorization.Solve(mat, vec);
      LinearEquation.Solve(mat, vec);

      return vec;
    }


  }
}
